#include "project.h"

double u1(double x);
double f1(double x);

int main(int argc, char* argv[])
{
    MultigridFactory M;
    const bool _VC_fl_registered = M.RegisterMultigrid({"VC", "full-weighting", "linear"}, CreateVC_fl);
    const bool _VC_fq_registered = M.RegisterMultigrid({"VC", "full-weighting", "quadratic"}, CreateVC_fq);
    const bool _VC_il_registered = M.RegisterMultigrid({"VC", "injection", "linear"}, CreateVC_il);
    const bool _VC_iq_registered = M.RegisterMultigrid({"VC", "injection", "quadratic"}, CreateVC_iq);
    const bool _FMG_fl_registered = M.RegisterMultigrid({"FMG", "full-weighting", "linear"}, CreateFMG_fl);
    const bool _FMG_fq_registered = M.RegisterMultigrid({"FMG", "full-weighting", "quadratic"}, CreateFMG_fq);
    const bool _FMG_il_registered = M.RegisterMultigrid({"FMG", "injection", "linear"}, CreateFMG_il);
    const bool _FMG_iq_registered = M.RegisterMultigrid({"FMG", "injection", "quadratic"}, CreateFMG_iq);
    
    std::tuple<std::string, std::string, std::string> schema;
    int n_grid, v1, v2, max_iteration;
    std::vector<double> boundary(2, 0.0);//boundary condition
    boundary[0] = u1(0);
    boundary[1] = u1(1);
    double w;//weighting coefficient
    double acc;//relative accuracy
    std::cin >> std::get<0>(schema) >> std::get<1>(schema) >> std::get<2>(schema) >>
	n_grid >> w >> v1 >> v2 >> max_iteration >> acc;
    Multigrid *m = M.CreateMultigrid(schema);
    m->set_n_grid(n_grid);
    m->set_RHS(boundary, f1);
    std::vector<double> u_precise(n_grid-1, 0.0);
    std::vector<double> u_compute(n_grid-1, 0.0);
    std::vector<double> error(n_grid-1, 0.0);
    std::vector<double> error_maxnorm(max_iteration+1, 0.0);
    std::vector<double> residual(n_grid-1, 0.0);
    std::vector<double> residual_L2norm(max_iteration+1, 0.0);
    for(int i = 0; i < n_grid-1; i++)
    	u_precise[i] = u1((i+1)*1.0/n_grid);
    double u_precise_maxnorm;
    u_precise_maxnorm = m->vector_max_norm(u_precise);
    for(int i = 0; i < n_grid-1; i++)
    	error[i] = u_precise[i] - u_compute[i];
    error_maxnorm[0] = m->vector_max_norm(error);
    for(int i = 0; i < n_grid-1; i++)
    	residual[i] = m->get_RHS()[i] - m->A_operater(u_compute)[i];
    residual_L2norm[0] = m->discrete_L2_norm(residual);
    double relative_accuracy = 0.0;
    relative_accuracy = error_maxnorm[0]/u_precise_maxnorm;
    int k = 0;
    if(std::get<0>(schema) == "VC")
    {
	std::vector<double> errorinit(n_grid-1, 0.0);
	while(k < max_iteration && relative_accuracy > acc)
	{
	    error = m->solver(errorinit, residual, w, v1, v2);
	    error_maxnorm[k+1] = m->vector_max_norm(error);
	    for(int i = 0; i < n_grid-1; i++)
		u_compute[i] = u_compute[i] + error[i];
	    for(int i = 0; i < n_grid-1; i++)
		residual[i] = m->get_RHS()[i] - m->A_operater(u_compute)[i];
	    residual_L2norm[k+1] = m->discrete_L2_norm(residual);
	    k = k + 1;
	    relative_accuracy = error_maxnorm[k]/u_precise_maxnorm;
	}
    }
    else
    {
	while(k < max_iteration && relative_accuracy > acc)
	{
	    error = m->solver(u_compute, residual, w, v1, v2);
	    error_maxnorm[k+1] = m->vector_max_norm(error);
	    for(int i = 0; i < n_grid-1; i++)
		u_compute[i] = u_compute[i] + error[i];
	    for(int i = 0; i < n_grid-1; i++)
		residual[i] = m->get_RHS()[i] - m->A_operater(u_compute)[i];
	    residual_L2norm[k+1] = m->discrete_L2_norm(residual);
	    k = k + 1;
	    relative_accuracy = error_maxnorm[k]/u_precise_maxnorm;
	}
    }
    std::vector<double> ratio_re(k, 0.0);
    for(int i = 0; i < k; i++)
	ratio_re[i] = residual_L2norm[i+1]/residual_L2norm[i];
    std::cout << "Multigrid method:" << std::get<0>(schema) << std::endl;
    std::cout << "Restriction operater:" << std::get<1>(schema) << std::endl;
    std::cout << "Interpolation operater:" << std::get<2>(schema) << std::endl;
    std::cout << "Parameter choice:" << std::endl;
    std::cout << "1.The number of grids in finest grid:" << n_grid << std::endl;
    std::cout << "2.The boundary condition:u0 = " << boundary[0] << ", u1 = " << boundary[1] << std::endl;
    std::cout << "3.The weighting coefficient of weighted Jacobi method:" << w << std::endl;
    std::cout << "4.The iterations before and v1 = " << v1 << ", v2 = " << v2 << std::endl;
    std::cout << "5.The max iterations:" << max_iteration << std::endl;
    std::cout << "6.The relative accuracy:" << acc << std::endl;
    std::cout << std::endl;
    std::cout << "residual_L2norm" << "  " << "reduction rate of the residual" << "  " << "error_maxnorm" << std::endl;
    for(int i = 0; i < k; i++)
	if(i == 0)
	    std::cout << residual_L2norm[i] << "                 " << "        " << "                    " << error_maxnorm[i] << std::endl;
	else
	    std::cout << residual_L2norm[i] << "                 " << ratio_re[i-1] << "                    " << error_maxnorm[i] << std::endl;
    return 0;
};

double u1(double x)
{
    double u;
    u = exp(sin(x));
    return u;
};

double f1(double x)
{
    double f;
    f = exp(sin(x))*cos(x)*(sin(x) - cos(x));
    return f;
};
